
#ifndef SHAPEANALYSIS_EXTRINSIC_DIRAC_H
#define SHAPEANALYSIS_EXTRINSIC_DIRAC_H

#include <vector>
#include <unordered_map>

#include <Eigen/Dense>
#include <Eigen/Geometry>

namespace gf {
    using std::vector;
    using std::unordered_map;
    using std::unordered_multimap;
    using Eigen::Quaternionf;
    using Eigen::Matrix;
    using Eigen::Vector;
    using face_id = unsigned int;
    using vertex_id = unsigned int;
    using edge_id = uint64_t;

    /**
     * @note this class is used to calculate dirac operator given a discrete polygon mesh
     * @note this algorithm may NOT be well performed given quadratic mesh or other mesh types,
     * @note i.e. mesh with faces which are not triangle.
     *
     * @cite for detailed information, please see [Ye 2018]https://onlinelibrary.wiley.com/doi/abs/10.1111/cgf.13494
     * @brief for usage, there will be a simple demo program in demo folder.
     * @brief in case you did not find it, it may NOT be finished, stay tuned!
     */
    class Extrinsic_Dirac {

    private:

        inline edge_id genIndex(const vertex_id idx_i, const vertex_id idx_j) {
            return (edge_id) ((edge_id) (idx_i << 32) + (edge_id) idx_j);
        }

        inline std::tuple<vertex_id, vertex_id> decomposeEdgeIndex(const edge_id e_idx) {
            return std::make_tuple((e_idx >> 32), (e_idx & 0xffffffff));
        }

    private:
        unsigned int num_vertex_per_face_;
        unsigned int num_faces_;
        unsigned int num_vertices_;
        unsigned int num_edges_;

        /// Willmore energy for mesh tesselation.
        float willmore_energy_;

        /// each edge correspond to a quaternion, in which real part is mean curvature of edge, imaginary part is vector of edge.
        /// in case two of adjacent faces are sharing a common plane,
        /// namely principle curvature will vanish, real part will equal to 0.
        vector<Quaternionf> mean_curvature_;

        /// areas for each face, for the purpose of Willmore energy calculation.
        vector<float> area_per_faces_;

        unordered_map<edge_id, unsigned int> mean_curvature_per_edge_;

        /// this map can be seen as a incident matrix of E2.
        /// incident matrix of E1 can be obtained using shift operators, thus no need to store them.
        unordered_multimap<edge_id, face_id> edge_faces_map_;
        unordered_map<face_id, unsigned int> area_per_face_;
        Matrix<Quaternionf, Eigen::Dynamic, Eigen::Dynamic> dirac_operator_;


    };
}


#endif //SHAPEANALYSIS_EXTRINSIC_DIRAC_H
